﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace InverseSqrtBenchmark {
    class Program {
        static void Main(string[] args) {
            const int Loops = 1000 * 1000 * 100;
            const float sq = 10.0f;
            Stopwatch watch = new Stopwatch();

            float isqrt = 0f;

            Console.WriteLine("findInverseSqrt(sq)");
            Console.WriteLine(Math.Abs(findInverseSqrt(sq) - 1.0f / (float) Math.Sqrt(sq)));
            watch.Reset();
            watch.Start();
            for (int i = 0; i < Loops; i++) {
                isqrt += findInverseSqrt(sq);
            }
            watch.Stop();
            Console.WriteLine(watch.Elapsed);

            Console.WriteLine("InvSqrt(sq)");
            Console.WriteLine(Math.Abs(InvSqrt(sq) - 1.0f / (float)Math.Sqrt(sq)));
            watch.Reset();
            watch.Start();
            for (int i = 0; i < Loops; i++) {
                isqrt += InvSqrt(sq);
            }
            watch.Stop();
            Console.WriteLine(watch.Elapsed);

            Console.WriteLine("Math.Sqrt(sq)");
            watch.Start();
            for (int i = 0; i < Loops; i++) {
                isqrt += 1.0f / (float) Math.Sqrt(sq);
            }
            watch.Stop();
            Console.WriteLine(watch.Elapsed);
        }
        public static float findInverseSqrt(float a) {
            const float eps = 1.0e-5f;
            float x = 2.0f / (a + 1);
            float ahalf = a * 0.5f;
            float residue;

            //while (Math.Abs(residue = x * (0.5f - ahalf * x * x)) > eps) {
            unsafe {
                int i = *(int*)&x;
                i &= 
                x += residue;
            }

            return x;
        }
        public static float InvSqrt(float x) {
            float xhalf = 0.5f * x;
            unsafe {
                int i = *(int*)&x; // store floating-point bits in integer
                i = 0x5f3759d5 - (i >> 1); // initial guess for Newton's method
                x = *(float*)&i; // convert new bits into float
            }
            x = x * (1.5f - xhalf * x * x); // One round of Newton's method
            return x;
        }
    }
}
